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ABSTRACT 

We report the results of new smoothed particle hydrodynamics (SPH) cal- 
culations of parabolic collisions between a subgiant or slightly evolved red-giant 
star and a neutron star (NS). Such collisions are likely to provide the dominant 
formation mechanism for ultracompact X-ray binaries (UCXBs) observed today 
in old globular clusters. In particular, we compute collisions of a 1.4 M NS 
with realistically modelled parent stars of initial masses 0.8 and 0.9 M Q , each at 
three different evolutionary stages (corresponding to three different core masses 
m c and radii R). The distance of closest approach for the initial orbit varies from 
r p = 0.04 R (nearly head-on) to r p = 1.3 R (grazing). These collisions lead to the 
formation of a tight binary, composed of the NS and the subgiant or red-giant 
core, embedded in an extremely diffuse common envelope (CE) typically of mass 
~ 0.1 to 0.3 Mq. Our calculations follow the binary for many hundreds of orbits, 
ensuring that the orbital parameters we determine at the end of the calculations 
are close to final. Some of the fluid initially in the giant's envelope, from 0.003 to 
0.023 M Q in the cases we considered, is left bound to the NS. The eccentricities of 
the resulting binaries range from about 0.2 for our most grazing collision to about 
0.9 for the nearly head-on cases. In almost all the cases we consider, gravitational 
radiation alone will cause sufficiently fast orbital decay to form a UCXB within 
a Hubble time, and often on a much shorter timescale. Our hydrodynamics code 
implements the recent SPH equations of motion derived with a variational ap- 
proach by Springel & Hernquist and by Monaghan. Numerical noise is reduced 
by enforcing an analytic constraint equation that relates the smoothing lengths 
and densities of SPH particles. We present tests of these new methods to help 
demonstrate their improved accuracy. 
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1. Introduction and Motivation 

Ultracompact X-ray binaries (UCXBs) are bright X-ray sources, with X-ray luminosities 
L > 10 36 erg s -1 , and very short orbital periods P < 1 hr. UCXBs are generally believed 
to be powered by accretion from a white dwarf (WD) donor onto a neutron star (NS), with 
the mass transfer driven by gravitational radiation (Rappaport et al. 1987). The abundance 
of bright X-ray binaries in Galactic globular clusters exceeds that in the field by many 
orders of magnitude, indicating that these binaries are formed through dynamical processes 
(Clark 1975). Indeed the stellar encounter rate in clusters has been shown to correlate 
with the number of close X-ray binaries (Pooley et al. 2003). In this paper, we use SPH 
calculations to investigate how UCXBs can be formed through direct physical collisions 
between a NS and a subgiant or small red giant (Verbunt 1987; Davies, Benz, & Hills 1992; 
Rasio & Shapiro 1991; Ivanova et al. 2005). In this scenario, the collision strips the subgiant 
or red giant, leaving its core orbiting the NS in an extremely diffuse common envelope 
(CE) of residual gas. Through the combined dissipative effects of CE evolution and, more 
importantly, gravitational radiation, the orbit decays until the core of the stripped giant 
overflows its Roche lobe and mass transfer onto the NS begins. 

Many recent studies of UCXBs have revealed their important role in a number of differ- 
ent contexts. They may be dominant in the bright end of the X-ray luminosity function in 
elliptical galaxies (Bildsten & Deloye 2004). UCXBs allow us to better understand in general 
the stellar structure and evolution of low-mass degenerate or quasi-degenerate objects (De- 
loye & Bildsten 2003; Deloye et al. 2005). They are also connected in a fundamental way to 
NS recycling and millisecond pulsar formation. Indeed, three out of five accretion-powered 
millisecond X-ray pulsars known in our Galaxy are UCXBs (Chakrabarty 2005). Finally, 
UCXBs may likely be the progenitors of the many eclipsing binary radio pulsars with very 
low-mass companions observed in globular clusters (Rasio et al. 2000; Freire 2005). 

Several possible dynamical formation processes for UCXBs have been discussed in the 
literature. Exchange interactions between NSs and primordial binaries provide a natural way 
of forming possible progenitors of UCXBs (Rasio et al. 2000). This may well dominate the 
formation rate when integrated over the entire dynamical history of an old globular cluster. 
However, it is unlikely to be significant for bright UCXBs observed today. This is because 
the progenitors must be intermediate-mass binaries, with the NS companion massive enough 
for the initial mass transfer to become dynamically unstable, leading to CE evolution and 
significant orbital decay. Instead, all MS stars remaining today in a globular cluster (with 
masses below the MS turn-off mass m to — 0.8 M & ) have masses low enough to lead to stable 
mass transfer (and the formation of wider LMXBs with non-degenerate donors). Alterna- 
tively, some binaries with stable mass transfer could evolve to ultra-short periods through 
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magnetic capture (Pylyser & Savonije 1988; Podsiadlowski et al. 2002). However, producing 
UCXBs through this type of evolution requires extremely efficient magnetic braking and a 
very careful tuning of initial orbital period and donor mass, and it is therefore very unlikely 
to explain most sources (van der Sluys et al. 2005a,b). 

Figure 1 helps provide the motivation for our consideration of collisions with subgiants or 
small red giants by showing that a significant fraction of collisions occur during these stages. 
This figure displays, for a star of mass O.9M , the stellar radius R as a function of time as 
well as the normalized number of physical collisions A(r) = r(t)dt/ r(t)dt, where to is 
the initial time of consideration, tf corresponds to the end of the red giant stage, r(t) is the 
rate of collisions, r(t) oc a(t) oc R 2 [1 + 2G(M + M NS )/(i?w^ )], a is the cross-section for the 
physical collision with a NS, M NS is the mass of the NS, and is the relative velocity at 
infinity. From the open circles in this figure we see, for example, that nearly 80% of the stars 
with M ~ O.9M that do collide will suffer their collision before leaving the main sequence 
(at t = 9.2 x 10 9 yr). This leaves more than 20% of such stars to collide while a subgiant 
or red giant. Of these, more than 60% collide while the star still has a radius R < 10R Q 
(see the triangular data point at the time when the radius R = 10R Q ). Therefore, due to 
their being larger than main sequence stars and evolving less rapidly than large red giants, 
subgiants and small red giants are therefore significant participants in stellar collisions. 

In a previous paper, we showed that direct physical collisions between NSs and subgiants 
or small red giants can provide a sufficient formation rate to explain the observed numbers 
of UCXBs (Ivanova et al. 2005). In this paper, we present in detail the hydrodynamics 
calculations of such collisions. Our paper is organized as follows. In §2, we present the SPH 
method. In §3, we present results from one- and three-dimensional test simulations that 
compare the performance of the variational and classical SPH equations of motion. In §4, 
we present the SPH models of the parent stars we consider, and in §5 these models are used 
in collision simulations. Finally in §6, we summarize our results and discuss future work. 

2. Numerical Methods 

SPH is the most widely used hydrodynamics scheme in the astrophysics community. It 
is a Lagrangian particle method, meaning that the fluid is represented by a finite number of 
fluid elements or "particles." Associated with each particle % are, for example, its position r^, 
velocity Vj, and mass mj. Each particle also carries a purely numerical smoothing length hi 
that determines the local spatial resolution and is used in the calculation of fluid properties 
such as acceleration and density. See Rasio & Lombardi (1999) for a review of SPH, especially 
in the context of stellar collisions. 
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Fig. 1. — Normalized number of physical collisions A shown for a star of mass O.9M in a 
cluster with = 10 km s -1 . The upper panel shows the evolution from the zero age main 
sequence to the top of the giant branch, and the lower panel zooms in to show the subgiant 
and red giant stages. The solid curve represents the stellar radius R (left axis). Symbols 
show A (right axis) at the current age: open circles are for a t corresponding to the zero 
age main sequence, while solid triangles are for a to corresponding to the terminal age main 
sequence. The markers (a), (b), and (c) in the bottom panel correspond to three of the stars 
modelled in this paper, namely SG0.9a, RG0.9b, and RG0.9c respectively (see §4). 
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Not surprisingly, there are many variations on the SPH theme. For example, one can 
choose to integrate an entropy-like variable, internal energy, or total energy. There are 
also different ways to estimate pressure gradient forces and a variety of popular choices for 
the artificial viscosity. Tests and comparisons of many of these schemes are presented in 
Lombardi et al. (1999). 

The classical formulations of SPH have proven more than adequate for numerous ap- 
plications, but because their underlying equations of motion implicitly assume that hi is 
constant in time, they do not simultaneously evolve energy and entropy strictly correctly 
when adaptive smoothing lengths are used (Rasio 1991; Hernquist 1993). In cosmological 
SPH simulations, for example, the resulting error in entropy evolution can significantly af- 
fect the final mass distribution (Serna, Dommguez-Tenreiro, & Saiz 2003). Seminal papers 
addressing the problem (Nelson & Papaloizou 1993, 1994; Serna et al. 1996) allowed for 
adaptive smoothing lengths, but in a somewhat awkward manner that was not generally 
adopted by the SPH community. 

More recently, Springel & Hernquist (2002) and Monaghan (2002) have used a varia- 
tional approach to derive the SPH equations of motion in a way that allows very naturally 
for variations in smoothing length. As we show in §3, this formalism works especially well 
when coupled with a new approach in which hi and pi are solved for simultaneously (Mon- 
aghan 2002). The idea is that some function of the particle density and smoothing length 
is kept constant (see §2.1), thereby satisfying the requirement of the variational derivation 
that the smoothing lengths be a differentiable function of particle positions. The purpose of 
this section is to present our implementation of these techniques. 



An estimate of the fluid density at rj is calculated from the masses and positions of 
neighboring particles as a local weighted average, 



2.1. Density and Smoothing Length 




(1) 
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where Wij(hi) = Wfcj, hi) is a smoothing (or interpolation) kernel with = |r s — r^j . We 
use the second-order accurate kernel of Monaghan & Lattanzio (1985), 

3 /r \ 2 3 /r \ 3 



W(r, /i) 



4 

r 

0, r > 2. 



Because the variational approach of deriving the SPH equations of motion requires that 
the smoothing lengths be a function of the particle positions, we choose an analytic function 
Hi(p) and solve hi = Hi(pi) simultaneously with equation (1); this idea is presented in 
Monaghan (2002) where it is credited to Bonet. Note that this solution can be found one 
particle at a time, as equation (1) depends only on the smoothing length hi and not on any 
hj (j 7^ i). The analytic function we use is 

H i (p) = {h^ i + b i p 1 /")- 1 , (3) 

where v is the number of dimensions, h maXj i is the maximum allowed smoothing length for 
particle i (the value approached as the density p tends to zero), and bi is a free parameter that 
can be adjusted according to the desired initial number of neighbors. For our one-dimensional 
test simulations in §3 that use equation (3), we set bi = 4/(mjA^v), where the constant 
Nn is a typical initial number of neighbors, and /i max ,i — +°° For those one- dimensional 
tests that do not implement equation (3), we simply adjust the smoothing lengths to always 
include the same number of neighbors N^, which is equivalent to enclosing the same neighbor 
mass as all particles are given the same mass. For our three-dimensional simulations, we 
choose h m3Xj i = 9600-R Q and bi to give approximately the desired initial number of neighbors, 
either N N w 32 (for N = 15,780) or N N w 48 (for N = 59,958); because /i maXji is much 
larger than the length scale of the problem, our results are not sensitive to its precise value. 
For computational efficiency, we solve for the smoothing lengths using a Newton-Raphson 
iterative scheme that reverts to bisection whenever an evaluation point would be outside the 
domain of a bracketed root. 



2.2. Pressure, Energy, and Entropy 

We associate with each particle i an internal energy per unit mass Ui in the fluid at iv 
In our one-dimensional test simulations, we implement a simple polytropic equation of state 
to determine pressure, 

Pi = (7 - , (4) 
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with the adiabatic index 7 = 5/3, corresponding to a monatomic ideal gas. We define the 
total entropy in the system as 



S = — !— - J~] rrii \n(uip\~ 
7 — 1 . ■*' 



(5) 



Although we need to consider only low-mass stars in our three-dimensional calculations 
of this paper, our code is quite general and includes a treatment of radiation pressure, which 
can be important in more massive objects: 



PikTj 1 4 

Pi = + o flT i > 

fa o 



(6) 



where k is the Boltzmann constant, a is the radiation constant, and fa is the mean molecular 
mass of particle i. The temperature T is determined by solving 



3 kTi alf 

Ui = - + 

2 fa pi 



(7) 



Equation (7) is just a quartic equation for Tj, which we solve via the analytic solution (see, 
e.g., Stillwell 1989) rather than through numerical root finding. In particular, suppressing 
the subscript i for convenience, the temperature is 

T= 1 -[-b+(b*-4cy/*}, 



where c = -q/(2b) +Y,b= (2F) 1 / 2 , Y = Y+ + Y_, 

Y 2\2 



Y ± = { q -± 
' ^ 16 



16 J \3aJ 



1/2' 



1/3 



and q = 3kp/(2afa). Careful programming of the expressions for T and Y_ can minimize 
roundoff error (see §5.6 of Press et al. 1992). 

Once the temperature T is known, the pressure of particle % is obtained from equa- 
tion (6). For our three-dimensional calculations, we define the total entropy in the system 
as 



E 



rrii 



3fc_ ln / 3^T -2/3 
2 fa \2 fa 1 



+ 



AaJf 
3pi 



(8) 
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2.3. Dynamic Equations and Gravity 

To evolve the system, particle positions are updated simply by = Vj and velocities 
by Vj = v| SPH ' ) + v l - Grav ' ) , where v| SPH ' 1 and v| Grav ' ) are the hydrodynamic and gravitational 
contributions to the acceleration, respectively. Various incarnations of the classical SPH 
equations of motion include 

„■ \Pi Pj J 



J 

and 



vf PH) = ~J2 m i (^iWijihi) + %V i W ij (h j ) + U ij V i W i3 \ . (10) 

j \Pi Pj ' 

The artificial viscosity (AV) term (see §2.4) ensures that correct jump conditions are 
satisfied across (smoothed) shock fronts. The symmetric weights Wij = Wji are calculated 
as n 

V«(hi) + w^)] . (ii) 



Springel & Hernquist (2002) derive a new acceleration equation for the case when an 
entropy-like variable is integrated. Monaghan (2002) generalized their work and showed 
that, even when integrating internal energy, the same acceleration equation should be used, 
namely 

vf PH) = (^o^iWM + V v . Wij(hj) + U .. Vi w ij ) , (12) 

j \ li iPi jPj J 

where 

* =i -it:2: m i-dh-- (13) 

The significant difference between equation (9) or (10) and equation (12) is the Q fac- 
tor. If the smoothing lengths do not change (Hi = hi = constant), then is unity and 
equation (12) reduces down to equation (10), one of the many classical variations of SPH. 
However, the deviation of Qj from unity corrects for errors that arise from a changing hi, 
and consequently both entropy and energy are evolved correctly. 

We use two MD-GRAPE2 boards to calculate the gravitational contribution to the 
particle acceleration as 

. (Grav) Grrij „ 

j^i V ^ 1 3 



where G is Newton's gravitational constant, is the softening parameter of particle i, and 
Tij is the unit vector that points from particle j toward particle i. For the SPH particles, q 
is set approximately to the initial smoothing length hi. After trying several different values 
of softening parameters for the core, we selected the one for each parent that yielded a gentle 
relaxation and accurate structure profiles (see Fig. 7). In all cases, the softening parameter 
of the core was comparable to but larger than the largest ti in the system. For the collision 
simulations, the softening parameter of the NS is then set to that of the core. 

The total gravitational potential energy is also calculated by the MD-GRAPE2 boards, 



which clearly excludes the gravitational self energy of the NS, core, and individual SPH 
particles. The GRAPE-based direct N-body treatment of gravity is crucial for maintaining 
excellent energy conservation even for very long runs. Using two MD-GRAPE2 boards, it 
takes about 2.5 seconds to calculate gravitational accelerations or potentials in our N = 
59, 958 simulations. 



Tests of various AV schemes are presented by Lombardi et al. (1999). For the one- 
dimensional shocktube tests of this paper, we implement the AV form proposed by Monaghan 
(1989) with a — (3 — 1 and rf = 0.01. For the collision calculations of this paper, we 
implement a slight variation on the form developed by Balsara (1995): 



as 




(15) 



2.4. Artificial Viscosity 




(16) 



where we use a — (3 — 1. Here, 



(Vi - Vj) • - Tj) fi + fj 



if ( Vi - Vj) ■ (Ti - Tj) < ; 



H>ij — 



ha(\ T i- T i\ 2 /hi3+v 2 ) 2 % 

o, 



if - w 3 



)-(r i -r j )>0, 



(17) 



where fi is the form function for particle % defined by 
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The factors if = 10 2 and rj' = 10 5 prevent numerical divergences. We calculate the 
divergence of the velocity field as 

v,) • ViWijihi) , (19) 

Vj) x ViWijihi) . (20) 

The function f\ acts as a switch, approaching unity in regions of strong compression (|V ■ 
v|j >> |V x v|j) and vanishing in regions of large vorticity (|V x v|j >> |V ■ v|j). Conse- 
quently this AV has the advantage that it is suppressed in shear layers. 

The only change in equation (16) from Balsara's original form is the inclusion of Q{ and 
flj, which equal one in simulations without adaptive smoothing lengths. Our experience is 
that the inclusion of these fl factors within allows for a more accurate AV scheme. 

2.5. Thermodynamics 

To complete the description of the fluid, Uj is evolved according to a discretized version 
of the first law of thermodynamics: 

7^ = E m i + \^) ^ - v i) • ViW«(^) . (21) 

We call equation (21) the "energy equation." A classical form of this equation has fij = 1 
and sometimes uses the symmetrized kernel Wij in place of Wij(hi). Although Monaghan 
(2002) uses for the AV term (only), we prefer equation (21) because it arises naturally 
in an introduction of AV for which p a /(Q a pl) is simply replaced by p a /(Q a pl) + n a &/2 in 
both the velocity and energy evolution equations (with a,b — i,j or We find that 

equation (21) treats shocks with essentially identical accuracy as the Monaghan formulation, 
and both formulations conserve total energy, X^ m i( v i ' + dui/dt) = 0. The derivation 
of equation (21) accounts for the variation of hi, so when we integrate it in the absence 
of shocks, the total entropy of the system is properly conserved even though the particle 
smoothing lengths vary in time. 

2.6. Integration in Time 

The evolution equations are integrated using a second-order explicit leap-frog scheme. 
For stability, the timestep must satisfy a Courant-like condition. Specifically, we calculate 



(V- v)i = ^-^rrijivj 



and the curl as 

(V x v)i = -E m i( v < 



Pi j 
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the timestep as 



At = Mini 



(Ai-i + A^)- 



For any SPH particle i, we use 

Afi,i = 

with 



C 



hi 



NA 



Max [Maxj (k^) , Maxj (/%)] 



K 



13 



Pi 



and 



Ai2,i — CiV,2 



in, 

2 ' 



ft 



1/2 



1/2 



(22) 



(23) 



(24) 



(25) 



| a i — ( a )j 

For the simulations presented in this paper, Cn,i = 0.6 to 0.9 and Cn,2 = 0.08 to 0.1. 
The Max,,- function in equation (23) refers to the maximum of the value of its expression 
for all SPH particles j that are neighbors with i. The denominator of equation (23) is an 
approximate upper limit to the signal propagation speed near particle i. The denominator 
inside the square root of equation (25) is the deviation of the acceleration of particle % from 
the local smoothed acceleration (a)j, given by 

m j a. j Wi j (hi) 



Pj 



(26) 



The advantage of including (a) i in this way is that the Lagrangian nature of SPH is preserved: 
the timestep would be unaffected by a constant shift in acceleration given to all particles. For 
the point particles % (namely, the core and the NS), we use At^ = CV^Min^ [r^/lvj — Vj|] 

1 /2 

and At2,« = CW^Miiij [r^/laj — a^]] , where the minimum is taken from all particles j 7^ i 
in the system. The incorporation of At2 enables us to use a larger Cn,i value and yields 
an overall more efficient use of computational resources. Indeed, we find this timestepping 
approach to be more efficient than any of those studied in Lombardi et al. (1999). 



2.7. Determination of the Bound Mass and Termination of the Calculation 

The iterative procedure used to determine the total amount of gravitationally bound 
mass Mi to the NS, M2 to the subgiant or red giant core, and M3 to the binary is similar to 
that described in Lombardi et al. (1996). As a minimal requirement to be considered part of 
a component, an SPH particle must have a negative total energy with respect to the center 
of mass of that component. More specifically, for a particle to be part of stellar component 
j = 1 (the NS) or j = 2 (the subgiant or red giant core), the following two conditions must 
hold: 
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1. vfj/2 + Ui — G(Mj — rrii)/dij must be negative, where is the velocity of particle i 
with respect to the center of mass of component j and d^ is the distance from particle 
i to that center of mass. 

2. dij must be less than the current separation of the centers of mass of components 1 
and 2. 

If conditions (1) and (2) above hold for both j = 1 and 2, then the particle is assigned to 
the component j that makes the quantity in condition (1) more negative. A particle % that 
is not assigned to component j = 1 or 2 is associated with component j = 3 (the common 
envelope) if only condition (1) is met for j = 1 or for j = 2. Remaining particles are assigned 
to the component 3 if they have a negative total energy with respect to the center of mass 
of the binary (mass Mi + M 2 ), or assigned to the ejecta otherwise. 

It should be noted that the method used in Ivanova et al. (2005) to determine mass 
components did not consider the possibility of a CE. Consequently the masses reported for 
the binary components were larger than here, although the conclusions of Ivanova et al. 
(2005) are unaffected. 

Once the mass components at a certain time have been identified, we calculate the eccen- 
tricity e and semimajor axis a of the binary from its orbital energy and angular momentum, 
under the approximation that the orbit is Keplerian. The kinetic contribution to the orbital 
energy comes simply from the total mass and momentum of each component in the center of 
mass frame of the binary. For the gravitational contribution, we first use the MD-GRAPE2 
boards to calculate the gravitational energy of all fluid in the union of components 1 or 2, 
subtract off the gravitational energy of just component 1, and then subtract the gravitational 
energy of just component 2. This way of calculating the orbital gravitational energy gives 
eccentricity e and semimajor axis a values that are much closer to being constant over an 
orbit than simpler methods that treated each star as a point mass. 

In all cases, we wait for at least 200 orbits before terminating a simulation. In many 
cases, in which orbital properties were still varying after 200 orbits or more than a few 
particles were still bound to the red giant core, we followed the evolution much longer (e.g., 
up to a total of 1743 orbits with over 1.2x 10 6 iterations for the case RG0.9b_RP3.82), in order 
to be certain that our orbital parameters were close to final. For comparison, computing 
power limited Rasio & Shapiro (1991) to simulate up to only 7 orbits in their SPH study of 
collisions between red giants and NSs. 



-13- 



3. Comparison of SPH Methods 

In this section, we compare the variational SPH scheme against the classical one for a 
few simple test cases. We perform both free expansion and shocktube tests in one dimension, 
problems that are particularly useful because of their known quasianalytic solutions. Indeed, 
free expansion of an initially uniform slab is just the limiting case of the Sod (1978) shocktube 
problem in which the density, pressure, and internal energy density in part of the tube are 
vansihingly small. In addition, we perform a free expansion test in three dimensions. 

To test the new dynamical equations (12) and (21), we implemented them in a one- 
dimensional SPH code and then compared their results for free expansion and shocktube 
problems with those of runs that use equations (12) and (21) but with = 1 for all particles. 
For each of these two formulations, we also vary the way in which the smoothing lengths 
are calculated: one run assigns hi such that each particle maintains the desired number of 
neighbors, and a second one implements the simultaneous solution of h-i and pi by using 
equation (3). Because we are integrating the internal energy equation, the level of overall 
energy conservation is limited by the integration technique and timestep, and not by the 
particular choice of equations being evolved: the increase in kinetic energy is always offset 
nicely by the decrease in internal energy. 

We begin by considering a simple case in which a uniform slab of 7 = 5/3 fluid, repre- 
sented by SPH particles possessing the same initial physical properties, is allowed to expand 
in one dimension into the surrounding vacuum. This fluid has an initial density p = 1 and 
pressure p = 1 over the range — 2 < x < 0. The initial specific energy values are set as 
Ui = pj ((7 — \)p) = 3/2, and the system is evolved with the help of the energy equation. In 
order to make a fair comparison of various methods, we turn off the AV and implement a 
constant timestep. 

We first examine how small the errors in energy and in entropy remain throughout the 
evolution. Figure 2 presents four calculations with N = 800 equal mass particles, N N = 5, 
and a timestep At = 10" 5 . From the symmetry of the curves in the bottom frame of Figure 2, 
we see that all four calculations yield an error 5(U + T) in total energy of essentially zero, 
which is not surprising given that an internal energy equation is being integrated. However 
the variational method with equation (3) (short dashed curve) does significantly better than 
any other method at keeping the errors in total internal energy U and total kinetic energy 
T individually close to zero. From the top frame, we see that the classical method (dotted 
curve) yields a spurious increase and then decrease in total entropy. When equation (3) is 
used in an otherwise classical method, the results are somewhat worse (dot-dashed curve). 
Similar results are obtained regardless of which exact form of the classical equations are 
used. Even the run that uses the new variational equations of motion but keeps the number 
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of neighbors strictly fixed (long dashed curve) produces unsatisfactory results: entropy tends 
to decrease with time. Merrily, implementing both the variational equations and equation 
(3) together (short dashed curve) produces the desired results of essentially zero net change 
in energy and entropy. Although the solution achieved by the other three approaches can 
be improved by choosing larger values of iV and Nn, we find that only the calculation that 
uses the variational equations and also simultaneously solves for hi and pi always has both 
excellent energy and entropy evolution. 

Figure 3 compares profiles of the specific internal energy u, density p, and velocity v 
against the analytic solution for two of the calculations treated in Figure 2: namely the 
classical method with fixed neighbor number (dot-dashed) and the variational method with 
simultaneously solution of particle densities and smoothing lengths (dashed). The classical 
calculation experiences significant numerical noise, manifested in the rapid oscillation of 
particles. Indeed, such oscilations are the root cause of the short timescale fluctuations seen 
in Figure 2 for the two calculations that do not simultaneously solve for hi and pi (the 
dotted and long-dashed curves). Although the variational method overshoots the analytic 
solutions for p and u near the edge of the rarefaction wave in Figure 3, it exhibits much 
smaller oscillations throughout most of the fluid and clearly represent a much more accurate 
solution. Interestingly, both solutions overestimate the density of the unperturbed fluid at 
— 1.8 < x < —0.2 by nearly 0.3%. This error, due to the relatively small number N N of 
kernel samplings for each pi calculation, is the same for all of these equally spaced particles 
that have not yet experienced either rarefaction wave. 

Energy conservation continues to be excellent in our tests even when the AV is turned on. 
For example when we use the same computational parameters as above with the variational 
equations of motion and the simultaneous solution of density and smoothing lengths, the 
total energy changes from its initial value of 3 to 2.999999998 by the time t = 0.15. As 
expected, energy conservation improves further when the timestep is decreased. The final 
error 5S in entropy over the same simulation is only 9 x 10~ 5 (as opposed to only 3 x 10~ 8 
when the AV is off). 

Given these encouraging results with free expansion tests, we now investigate how the 
variational equations and the simultaneous solution of hi and pi work for non-adiabatic 
processes. In the standard Sod (1978) shocktube test, two fluids, each with different but 
uniform initial densities and pressures are placed next to each other and left to interact 
over time. The initial conditions for our shocktube tests are similar to those in Rasio & 
Shapiro (1991), except that here the adiabatic index 7 = 5/3. In particular, one slab of 
fluid has a denisty pi — 1 and pressure pi — 1 over the range — 1 < x < 0, while a second 
slab has denisty p r = 0.25 and pressure p r = 0.12402 over the range < x < 1. From 
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Fig. 2. — As a function of time t, the errors in total entropy S, total internal energy U, and 
total kinetic energy T for one dimensional free expansion calculations with N = 800 and 
Nn = 5. Each integration method is represented by a different curve: classical with fixed 
neighbor number (dotted), classical with simultaneous solving of hi and pi (dot-dashed), 
variational with fixed neighbor number (long dashed), and variational with simultaneous 
solving of hi and pi (short dashed). In the bottom panel, the curves that are negatively 
valued at early times are for ST, while the others are SU. 
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Fig. 3. — The left column shows the specific internal energy u, density p, and velocity v 
plotted as functions of position x at time t — 0.15 for two of the calculations featured in Fig. 
2: purely classical (dotted) and variational with equation (3) (dashed). The solid curve in 
the left column gives the quasianalytic solution, while the right column gives the deviations 
from that solution. 
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the quasianalytic solution for these initial conditions, the entropy should increase at a rate 
dS/dt = 0.03834, while the energies change at a rate dT/dt = —dU/dt = 0.9751, where we 
have included the effects both of the shock and of the free expansion at the edges. 

The energy equation is integrated, with initial values of specific internal energy set, 
except for the particles on the edges of the system, according to Ui = PPi~ 1 p~^/('J — 1), 
where pi is calculated from equation (1), p and p are either pi and pi or p r and p r (depending 
on whether the particle is at a negative or positive x position), and ( = log(p r /pi) / \og(p r / pi) 
is chosen so that Ui varies smoothly across x = 0. The particles near an edge (near x — ±1) 
are given the same Ui as the bulk of the particles in that slab of fluid. As a result, Ui is 
constant on each side, except over a transition region of a few smoothing lengths near x — 0. 

Our shocktube calculations implement the same four formulations as in the free expan- 
sion tests. We employ N = 800 equal mass particles with = 5. To help make a fair 
comparison of methods, we set a constant timestep At = 10~ 5 . As with the free expansion 
tests, we begin by examining the error in the energies and entropy as a function of time. 
From the symmetry in the curves of the bottom panel of Figure 4, we see that ST = —5U 
to a high level of precision, as expected. That is, all four cases yield an error 5(U + T) in 
total energy of essentially zero. We note that implementing both the variational equations 
and equation (3) together (short dashed curve) produces the smallest errors in T and U, 
a result that continues to hold even for other choices of N, Nn, and /i max ,j- From the top 
two panels of this figure, we see that the entropy evolution is also most accurate when the 
variational equations of motion are used and the particle densities and smoothing lengths 
are found simultaneously. The error in entropy that does remain could be further reduced 
by a more sophisticated AV scheme. The next most accurate solution uses the variational 
equations with fixed neighbor number. Neither of the two classical formulations yield an 
accurate solution for these computational parameters. 

Figure 5 compares profiles against the analytic solution for the cases of the purely 
classical method (dotted curve) and the variational with equation (3) method (dashed curve). 
The purely classical method contains considerable noise, especially in the rarefaction (—0.2 < 
x < —0.05). These fluctuations can be diminished with a stronger AV, but at the expense of 
further inaccuracy in the total entropy in the system. In contrast, the dashed curve displays 
very little noise or error throughout most of the profiles, which is a direct result of solving 
for particle density and smoothing length simultaneously. Indeed, regardless of the choice 
of N, -/Vjv, and h maXti , the variational formulation with the simultaneous solution of hi and 
Pi yields results that are as good as, or often significantly better than, that from any other 
formulation considered. 

Figure 6 helps show that the SPH method in three dimensional calculations can also 
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Fig. 4. — Change in entropy AS, error in entropy SS, error in kinetic energy ST, and error 
in internal energy SU for a set of one dimensional shocktube test calculations with N = 800 
and = 5. As in Fig. 2, each integration method is represented by a different curve: 
classical (dotted), classical with simultaneous solution of hi and pi (dot-dashed), variational 
(long dashed), and variational with simultaneous solution of hi and pi (short dashed). The 
solid curve in the top panel represents the quasianalytic solution. In the bottom panel, the 
curves that are positive at early times are SU, while the remaining curves are ST. 
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Fig. 5. — Profiles at a time t — 0.15 for two of the one dimensional shocktube test calcula- 
tions featured in Fig. 4: classical with fixed neighbor number (dotted) and variational with 
simultaneous solution of hi and p i (dashed). 
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be improved by the combination of the variational equations of motion and the idea of 
constraining p { and hi analytically. In this test, we turn off AV, set vj Grav ^ = for all particles 
i, use a constant timestep of about 0.02 hours, and consider the free expansion of a spherically 
symmetric gas distribution, namely the RG0.9b red giant model presented in §4 and used in 
some of the collision simulations of §5. When the classical acceleration and internal energy 
evolution equations are integrated, the total entropy of the system spuriously increases with 
time. Figure 6 demonstrates this behavior for the case where f2j = 1 in equations (12) and 
(21). In contrast, the variational equations allow the entropy to be properly conserved (see 
the dashed line in the top panel). 



4. Modeling the Parent Stars 

In this section we present our procedure for modelling the parent stars that are used 
in the collision simulations of §5. Because the size of the NS is many orders of magnitude 
below our hydrodynamic resolution, we take the usual approach of modeling it as a point 
mass that interacts gravitationally, but not hydrodynamically, with the rest of the system. 
For our (sub)giant models, we first use a stellar evolution code developed by Kippenhahn et 
al. (1967) and updated as described in Podsiadlowski et al. (2002) to evolve stars of mass 
M = 0.8 and O.9M , primordial helium abundance Y = 0.25, and metallicity Z = 0.001, 
without mass loss. These stellar models were computed using OPAL opacities (Rogers & 
Iglesias 1992) and supplemented with opacities at low temperatures (Alexander & Ferguson 
1994). For each mass, we considered models corresponding to different evolution stages and 
hence different core masses. Table 1 gives the initial properties of the seven O.8M and 
0.9M© subgiant and red giant SPH models, including one higher resolution parent model. 
Column (1) gives the name of the parent star model, while Column (2) gives its mass. The 
next two columns give data resulting from the stellar evolution calculation: column (3) gives 
the stellar radius R of the parent, and column (4) lists the core mass. The last two columns 
present parameters, discussed more below, that are relevant to the SPH realization of each 
model: namely, column (5) shows the mass of the central gravitational point particle and 
column (6) lists the spacing of the hexagonal close-packed (hep) lattice cells. The models 
we are considering are either subgiants (i.e., in the Hertzsprung gap) or small red giants 
(i.e., near but after the base of the red giant branch). As discussed in Ivanova et al. (2005) 
and shown in Figure 1, such stars are the most likely to collide. In order to generate SPH 
models of these subgiant and red giant stars, we initiate a relaxation run by placing SPH 
particles on an hep lattice out to a radius of approximately two smoothing lengths less 
than the stellar radius calculated by the evolution code. Particle masses are assigned to 
yield the desired density profile. Although in general the use of unequal mass particles can 
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Fig. 6. — Entropy S (in cgs units), kinetic energy T, internal energy U, and total energy 
T + U in the three dimensional free expansion of a O.9M red giant star after gravity is 
turned off. The dotted curve results from the classical SPH equations, while the dashed 
curve is from the variational equations. In the bottom panel, the T curves increase with 
time, the U curves decrease, and the total energy T + U is essentially constant regardless of 
the method. In both calculations, equation (3) is used to set the smoothing lengths. 
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Table 1: Parent Star Characteristics 



Parent Star 


M 


R 


/n 
' ' v c 


777. 








(Re) 


\ kJ J 




(Re) 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


SG0.8a 


0.8 


1.60 


0.10 


0.15 


0.10 


RG0.8b 


0.8 


3.17 


0.19 


0.24 


0.20 


RG0.8c 


0.8 


4.43 


0.22 


0.27 


0.28 


RG0.8c_hr 


0.8 


4.43 


0.22 


0.25 


0.19 


SG0.9a 


0.9 


2.02 


0.12 


0.18 


0.13 


RG0.9b 


0.9 


5.31 


0.23 


0.28 


0.33 


RG0.9c 


0.9 


6.76 


0.25 


0.29 


0.43 



lead to spurious mixing in SPH calculations (e.g., Lombardi et al. 1999), we do not expect 
such effects to be significant in the collision simulations of this paper. Because the NS is 
modelled as a pure point particle, there can be no mixing of the SPH particles between the 
two stars. Furthermore, the entropy gradients in the (sub)giant parent stars and especially 
in the post-collisional shocked fluid will suppress spurious mixing. 

We choose an hep lattice for its stability (Lombardi et al. 1999). Each primitive hexag- 
onal cell has two lattice points and is spanned by the vectors ai, a 2 , and a 3 , where &i and a 2 
have an included angle of 120° between them and magnitudes a\ — a 2 , and a 3 is perpendic- 
ular to ai and a 2 with magnitude a 3 = 1.633 ai (see, e.g., Kittel 1986). In six of the parent 
models, we use N = 15, 780 SPH particles and ~ 32 to model the gaseous envelope of 
the parent. For the the most evolved O.8M parent, we also generated a model, RG0.8c_hr, 
using N = 59, 958 SPH particles and N N 48 neighbors, in order to help estimate the 
uncertainty in our final orbital parameters. 

Because the density in the core of these parent stars is roughly 10 5 to 10 s times larger 
than their average density, it is not possible to model the cores with SPH particles. Instead, 
we model the core as a point mass that interacts gravitationally, but not hydrodynamically, 
as suggested by Rasio & Shapiro (1991), among others. The core point mass is positioned 
at the origin with the six SPH particles nearest to it located on lattice points at a distance 
of ai/^/2. We set the mass of the core point mass m pt by subtracting the sum of the SPH 
particle masses from the known total mass of the star. For our parent models, m pt ranges 
from 0.15M© (SG0.8a) to O.29M (RG0.9c). Because the resolution in our parent models is 
much larger than the size of the core, the core point mass m p t is always somewhat larger than 
the physical core mass m c , usually by about 0.05M Q for our iV = 15, 780 models (see Table 
1). The increased resolution of model RG0.8c_hr allows for the point mass m pt = 0.25M Q to 
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be closer to the actual core mass m c = O.22M given by our stellar evolution code than it is 
in the RG0.8c model. 

After the initial parameters of the particles have been assigned, we relax the SPH fluid 
into hydrostatic equilibrium. During this process, we adjust the smoothing lengths and bi 
values at each iteration in order for each particle to maintain approximately the desired 
number of nearest neighbors. We employ both artificial viscosity and a drag force to assist 
with the relaxation. We also hold the position of the core point mass fixed. Figure 8 plots 
SPH particle data for one of our relaxed parent models, RG0.9c. Although the core of the 
star cannot be resolved in the innermost ~ 2 smoothing lengths, the thermodynamic profiles 
of the SPH model nicely reproduce those from the stellar evolution code. 



5. Results of Collision Simulations 

In this section, we report on the results of 32 simulations of parabolic collisions between 
a subgiant or small red giant and a 1.4M NS. The results of these collisions are summarized 
in Table 2. In terms of computational time, most of the N = 15, 780 particle runs lasted 
somewhere in the range from 5 to 20 days when run on a typical workstation with a Pentium 
III or IV processor. As for the three N = 59, 958 particle runs, the r p = O.96i? case ran 
for four weeks and completed its 204 orbits in nearly 7.9 x 10 4 iterations, the r p = 1.91i? 
case lasted nearly 6 weeks and completed 219 orbits in more than 1.2 x 10 5 iterations, and 
the r p = 3.82_R case took 10 weeks to complete 450 orbits in about 3.25 x 10 5 iterations. 
When multiple simulations were running, the time per iteration would vary depending on 
the availability of the MD-GRAPE2 boards: the hydrodynamics was done on separate work- 
stations, but communication with the same host of the MD-GRAPE2 boards was necessary 
for all calculations. 

Figure 9 presents particle plots for one of these collisions, case RG0.9b_RP3.82, in which 
a O.9M red giant collides with a 1.4M NS at a periastron separation of 3.82i? and an initial 
separation of 48i? . Frame (a) shows the two bodies just prior to impact. Frame (b) takes 
place during the first pericenter passage, while frame (c) shows the first apocenter passage. 
At this time, much of the mass originally bound to the red giant star has been transferred to 
the NS. Frame (d) occurs during the second pericenter passage, and, approximately twenty 
hours later, frame (e) occurs during the second apocenter passage. Finally, frame (f) shows 
a snapshot from the sixth apocenter passage, at a time of 277.5 hr. At this late time, 
few particles remain bound to the subgiant core. During all of the stellar collisions we 
considered, some SPH particles form a CE that surrounds the binary. Various projections 
of the SPH particles, including those particles in the CE, can be seen in Figure 10 at two 
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Fig. 7. — Pressure p, density p, temperature T, and mean molecular mass p versus radius r 
for the parent stars with initial masses of O.8M (left column) and O.9M (right column), 
as determined by our stellar evolution code. The three curves in each plot represent the 
different parent stars, from least evolved (solid line) to most evolved (dotted line). 
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Fig. 8. — Properties of a 0.9M© parent star model (RG0.9c) as a function of radius. In 
the left column (from bottom): mean molecular mass, temperature, density, and base 10 
logarithm of pressure. The dashed curves in these four graphs represent the calculated profile 
of the star, while the points represent our SPH particles. In the right column (from bottom): 
radial gravitational acceleration (lower curve) and hydrodynamic acceleration (upper curve), 
smoothing length, and individual particle mass. 
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Fig. 9. — Two-dimensional projections onto the orbital plane of SPH particle and point- 
mass positions at various times in case RG0.9b_RP3.82. The red giant core is illustrated 
with a solid circle and the NS with a solid triangle. Particles are colored according to the 
component to which they are currently considered bound: blue for the NS, red for the red 
giant core, green for the CE, and magenta for the ejecta. 
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different times for case RG0.9b_RP3.82. At 277.5 hr, there are 158 SPH particles considered 
bound to the NS, 3 SPH particles bound to the subgiant core component, and 4631 particles 
in the CE. By a time of 12,380 hr, when the simulation is terminated, these numbers have 
stabilized to 1, 0, and approximately 600, respectively. Although the outer layers of the CE 
are expanding even at the end of the simulation, this fluid is still gravitationally bound to 
the binary, with a mass of approximately 0.06M Q in this case. 

Figures 11 and 12 show energies as a function of time for the cases RG0.9b_RP3.82 and 
RG0.8c_RP0.96, respectively. The sharp oscillating peaks in the total kinetic and potential 
energy plots correspond to pericenter passages, which become more frequent as the collision 
progresses. During the late times shown in the right column, the binary has stabilized, and 
thus the internal energy and frequency of kinetic and potential energy oscillations remain 
relatively constant. 

For our N = 15, 780 runs, the difference between the maximum and minimum total 
energy is typically ~ 10 46 erg, at most a few percent of the total energy in the system. Energy 
conservation in our N = 59, 958 runs is better by a factor of ~ 2. We find that decreasing 
Cjv,i and Cjv,2 makes energy be conserved even more accurately, but does not significantly 
affect final masses or orbital parameters. Energy conservation at late times, when the binary 
is left surrounded by a CE, is typically excellent. The level of energy conservation tends to 
improve somewhat as we consider more evolved parent stars or larger impact parameters. In 
all of our runs, angular momentum conservation holds at an extremely high level of accuracy, 
typically at the 10~ 4 % level or better. Even in the case of RG0.8c_RP5.73, one of our longest 
runs, angular momentum is conserved to better than 3 x 10~ 4 %. 

Figure 13 shows component masses and orbital parameters as a function of time for case 
RG0.9b_RP3.82. By 200 hours (during the third pericenter passage), the subgiant (M_) has 
been nearly completely stripped. This lost mass is transferred to the NS (Mi), transferred 
to the CE (M 3 ), or, most commonly, ejected to infinity. By the end of the collision, Mi is 
approximately 1.4O3M and the eccentricity e has stabilized near 0.43. Figure 14 shows the 
evolution of the component masses and orbital parameters for the collision RG0.8b_RP3.82, 
a case that we carried out to over 1000 orbits. The perturbation to the orbit near a time 
of t = 1700 hours, most easily seen in the eccentricity plot, occurs when the final gas is 
stripped from the subgiant core. After this time, the orbital parameters for the inner binary 
are essentially unaffected, although there is a very gradual decrease in the mass M 3 considered 
to be in the CE. Figure 15 shows the same parameters for the collisions RG0.8c_RP0.96 and 
RG0.8c_hr_RP0.96. These two cases both involve the same, evolved parent star, but the 
latter case uses a model with more SPH particles (N = 59, 958 as opposed to AT — 15, 780) 
and a less massive central point particle (m pt = O.25M as opposed to m pt = O.27M ). In 
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Fig. 10. — Two dimensional projections of SPH particle and point-mass positions at two 
select times in the case RG0.9b_RP3.82: the left column is at 277.5 hr, as in frame (f) of 
Fig. 9, and the right column at 12,380 hr. Particle colors are as in Fig. 9. 
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Fig. 11. — Energies as a function of time for collision RG0.9b_RP3.82. The left column 
presents the time evolution up to 700 hr. The right column presents a 60 hr time interval 
much later in the simulation. The gravitational potential energy W curve is the most nega- 
tive, the total energy E curve is nearly horizontal, the internal energy U curve approaches a 
small positive energy at late times, and the kinetic energy K curve is the positive one with 
large fluctuations that are synchronized with the phase of the orbit. These energies do not 
include the self-energy of the core or NS. 
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Fig. 12.— Like Figure 11, but for the collision RG0.8c_RP0.96. 
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Fig. 13. — Masses Mj, semimajor axis a, and eccentricity e as a function of time for collision 
RG0.9b_RP3.82. The masses M 1 , M 2 , and M 3 correspond to the masses containing NS, the 
subgiant core, and the CE gravitationally bound to the binary system, respectively. The 
dashed vertical lines correspond to the moments of 25, 50, 100, 200, 400, 800, and 1600 
completed orbits. 
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Fig. 14. — Like Fig. 13, but for the collision RG0.8b_RP3.82. The dashed vertical lines 
correspond to the moments of 25, 50, 100, 200, 400, and 800 completed orbits. 
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both cases, the relatively small periastron separation leads to an eccentric orbit in which 
the red giant (M 2 ) is stripped of most of the initial mass within 50 hours, leaving only 
the point particle. As the binary stabilizes, Mi remains constant at 1.4O3M in both cases. 
Although the details of the M 3 curves differ at early and intermediate times, both cases settle 
to essentially the same final common envelope mass M3 of about 0.08 or O.O9M . In the 
case RG0.8c_RP0.96, the eccentricity e tends toward 0.71 and the semimajor axis a stabilizes 
near 1.34R at late times, while in the case RG0.8c_hr_RP0.96, these quantities aproach 0.62 
and 1.15i? , respectively. These differences in final orbital parameters are largely due to the 
different m pt values. For a given total mass and impact parameter, the trend is for smaller m pt 
to yield somewhat less eccentric and tighter orbits: this is true regardless of whether the m pt 
is smaller because of increased numerical resolution (compare the higher resolution collisions 
of RG0.8c_hr models with their medium resolution counterparts) or because the parent is in 
a different evolutionary stage (compare, for example, RG0.8c_RP0.96 and SG0.8a_RP0.96). 

By considering plots like those of Figures 13, 14, and 15, we are able to determine final 
component masses and orbital parameters for all of our calculations (see Table 2 and Fig. 
16). 

6. Discussion and Future Work 

We have shown that collisions between NSs and subgiants or small red giants naturally 
produce tight eccentric NS-WD binaries that can become UCXBs. These results reinforce 
our findings in Ivanova et al. (2005), where we consider the rate of UCXB formation and 
show that all currently observed UCXBs could be explained from this channel. In addition, 
our paper has demonstrated the importance in SPH codes of using the variational equations 
of motion and simultaneously solving for particle densities and smoothing lengths. 

We have found no significant complications applying the variational equations of motion 
or the simultaneous solution of pi and hi in either our one- or three-dimensional hydrody- 
namics code. Our simple tests of these methods in §3 help demonstrate the significantly 
improved accuracy of SPH when these approaches are applied together. In particular, the 
combination of these approaches produces a change in entropy that most closely matches 
that of known solutions, without disturbing energy conservation. Additional shocktube tests 
show similar results. The reduction in numerical noise that occurs when particle densities 
and smoothing lengths are found simultaneously is due to the resulting gradual variation of 
smoothing lengths from iteration to iteration. Note that if the smoothing lengths are instead 
determined by requiring that a particle have some desired number or mass of neighbors, then 
these smoothing lengths can typically take on a range of values. Noise is then introduced by 
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Fig. 15. — Like Fig. 13 and 14, but for two similar collisions: the thin curve is for 
RG0.8c_RP0.96 (N = 15,780 SPH particles and m pt = O.27M ), while the bold curve is 
for RG0.8c_hr_RP0.96 (N = 59,958 and m pt = O.25M ). 
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Fig. 16. — Mass M 3 of the CE, mass Mi of the NS and its bound fluid, final semimajor axis 
a, and final eccentricity e as a function of normalized periastron separation r p /R. The left 
column (open shapes) corresponds to collisions with a O.8M subgiant or giant star and the 
right column (solid shapes) corresponds to collisions with a O.9M subgiant or giant star. 
From least to most evolved: SG0.8a and SG0.9a are represented with triangles, RG0.8b and 
RG0.9b with squares, and RG0.8c and RG0.9c with circles. The symbol ® refers to collisions 
involving our RG0.8c_hr model. 
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Fig. 17. — Dependence of the gravitational radiation merger time r gw on post-collision semi- 
major axis a and eccentricity e. The points types are as in Fig. 16. The area of the shape 
is indicative of r/(dR/dt), where r is the rate of collisions for a star of radius R, and hence 
is larger for the less evolved subgiants, which are more likely to be involved in a collision. 
The dashed curves are curves of constant gravitational merger time, as labelled, for a 1.4 
M Q NS and a O.25M WD. The orbital period P values on the right axis are for these same 
masses and are calculated simply from Kepler's third law. The hatched area shows how 
the r gw = 10 10 yr curve changes when we vary the binary parameters: the upper boundary 
corresponds to a 1.5 M Q NS with a 0.45 M Q WD, and the lower boundary corresponds to a 
1.3 M Q NS with a 0.15 M WD. 
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abrupt changes in the smoothing length as neighbors enter or leave the kernel. 

In our parabolic collisions of a NS and a subgiant or red giant, we find the stripping of 
the core to be extremely efficient. In all but two cases, which were both grazing collisions 
with an evolved O.8M red giant, the core was stripped completely of SPH fluid. We could 
not, however, resolve the innermost fluid outside the core even in our parent models, and so 
the best we can do is place an upper limit of m pt — m c (usually ~ O.O5M ) on the actual 
amount of mass that would remain bound to the core after the collision. As the radius of the 
core is about O.O2i? for all of our (sub)giant parent stars, to resolve the fluid down to the 
core would require an increase in the number of particles by roughly a factor (O.O2i? /ai) 3 , 
where the a\ values that set the particle spacing are given in the last column of Table 1. 
Therefore to get the core mass correct for our least evolved O.8M star would require about 
2 x 10 6 particles. The other parent stars would require more particles, up to about 2 x 10 s 
for our most evolved O.9M RG. 

Some of the fluid initially in the sub giant or red giant envelope, from 0.003 to O.O23M 
in the cases we considered, is left bound to the NS. This fluid will likely be accreted onto 
the NS, potentially recyling it to millisecond periods. By carrying out our collisions between 
a NS and a subgiant or red giant to many orbits, we were able to identify the existence 
of a residual, distended CE surrounding a binary consisting of the NS and the subgiant or 
giant core. The ultimate fate of this diffuse CE is rather uncertain. This gas will likely 
be quickly ejected by the radiation released by accretion onto the NS. Nevertheless, even 
a brief CE episode phase will only increase the rate of orbital decay as compared to that 
from gravitational radiation alone. The gravitational merger timescale therefore provides an 
upper limit on how long it will take for a UCXB to form. 

When we apply the Peters (1964) equations to these post-collision systems, we find that 
most of them inspiral on rather short timescales (Fig. 17). We conclude that all collisions 
between a subgiant and a NS, as well as all but the most grazing RG-NS collisions, can 
produce UCXBs within a Hubble time. As seen in Figure 17, high eccentricities are an 
important factor in keeping merger times short. However, as discussed in Ivanova et al. 
(2005), even if all binaries were able to circularize quickly (compared to the gravitational 
merger time), a large fraction of post-collision systems would still merge in less than the 
cluster age, as can be seen directly from Figure 17. 

If we could run all of our simulations at much higher resolution, we would expect the 
details of the resulting orbital parameters to change but our conclusions to be unaffected. In 
particular, higher resolution calculations would allow the point mass m pt to be closer to the 
true core mass m c and could therefore more accurately determine the mass left bound to the 
(sub)giant core. Using our three high resolution calculations as a guide, we expect that the 
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resulting eccentricity e and semimajor axis a values of the binary would also be somewhat 
decreased. However, judging from the shape of the constant r gw curves in Figure 17 and the 
change in e and a that resulted by going from N = 15, 780 to N = 59, 958, we would not 
expect the gravitational merger time to be dramatically affected. That is, we would still 
typically be left with a tight, eccentric binary that would merge within a Hubble time. 

Bildsten & Deloye (2004) have recently shown that the cutoff in the observed luminosity 
function of extragalactic LMXBs can be explained if nearly all of the UCXB progenitors 
consisted of a NS accreting from a compact object with a mass comparable to the helium 
core mass at turnoff. Our scenario for UCXB formation provides a natural explanation for 
such progenitors. When our channel for UCXB formation is followed, the donor will be a 
helium WD, perhaps with small amounts of hydrogen left over from the (sub)giant envelope. 
It has long been known that the 11.4 minute UCXB 4U 1830-20 contains exactly such a 
helium donor (Rappaport et al. 1987). Recently, in 't Zand et al. (2005) have also identified 
2S 0918-549 as a likely UCXB with a helium WD donor. Although a subset of UCXBs may 
contain C-0 WD donors (see Juett & Chakrabarty (2005) and references therein), these 
systems could result from collisions between a NS and an assymptotic giant branch star. 

We thank R. Bi, S. Fleming, E. Gaburov, M. Rosenfeld, and the anonymous referee for 
helpful comments and contributions. This material is based upon work supported by the 
National Science Foundation under Grants No. 0205991 and 0507561. N.I. and F.A.R. also 
acknowledge support from NASA Grants NAG5-12044 and NNG04G176G at Northwestern 
University. This work was also partially supported by the National Computational Science 
Alliance under Grant AST980014N. 
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Table 2: Summary of Collisions 



Collision 


r p 


AM C 


AM NS 


M 3 


e 


a 


-^orbit 




(Re) 




(Mq) 






(Rq) 




(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


SG0.8a_RP0.24 


0.24 


0.0 


0.013 


0.23 


0.85 


0.24 


525 


SG0.8a_RP0.48 


0.48 


0.0 


0.013 


0.20 


0.80 


0.28 


550 


SG0.8a_RP0.96 


0.96 


0.0 


0.013 


0.24 


0.64 


0.33 


334 


RG0.8b_RP0.24 


0.24 


0.0 


0.005 


0.15 


0.89 


0.65 


302 


RG0.8b_RP0.48 


0.48 


0.0 


0.005 


0.11 


0.82 


0.77 


235 


RG0.8b_RP0.96 


0.96 


0.0 


0.005 


0.11 


0.70 


0.98 


207 


RG0.8b_RP1.91 


1.91 


0.0 


0.005 


0.10 


0.54 


1.26 


259 


RG0.8b_RP3.82 


3.82 


0.0 


0.023 


0.07 


0.33 


2.28 


1038 


RG0.8c_RP0.24 


0.24 


0.0 


0.003 


0.15 


0.89 


0.95 


201 


RG0.8c_RP0.96 


0.96 


0.0 


0.003 


0.09 


0.71 


1.34 


281 


RG0.8c_RP1.91 


1.91 


0.0 


0.003 


0.10 


0.64 


1.78 


298 


RG0.8c_RP3.82 


3.82 


0.0 


0.003 


0.06 


0.34 


2.19 


862 


RG0.8c_RP5.73 


5.73 


0.003 


0.010 


0.05 


0.18 


2.70 


1351 


RG0.8c_hr_RP0.96 


0.96 


0.0 


0.003 


0.08 


0.62 


1.15 


204 


RG0.8c_hr_RP1.91 


1.91 


0.0 


0.005 


0.09 


0.51 


1.62 


219 


RG0.8c_hr_RP3.82 


3.82 


0.003 


0.005 


0.08 


0.18 


2.18 


450 


SG0.9a_RP0.24 


0.24 


0.0 


0.013 


0.28 


0.92 


0.29 


371 


SG0.9a_RP0.48 


0.48 


0.0 


0.013 


0.22 


0.88 


0.33 


295 


SG0.9a_RP0.96 


0.96 


0.0 


0.013 


0.34 


0.64 


0.44 


539 


SG0.9a_RP1.43 


1.43 


0.0 


0.013 


0.22 


0.53 


0.41 


288 


SG0.9a_RP1.67 


1.67 


0.0 


0.013 


0.29 


0.54 


0.46 


330 


RG0.9b_RP0.24 


0.24 


0.0 


0.003 


0.18 


0.91 


0.97 


297 


RG0.9b_RP0.48 


0.48 


0.0 


0.003 


0.11 


0.85 


1.07 


531 


RG0.9b_RP0.96 


0.96 


0.0 


0.003 


0.14 


0.72 


1.26 


365 


RG0.9b_RP1.91 


1.91 


0.0 


0.003 


0.12 


0.55 


1.65 


434 


RG0.9b_RP3.82 


3.82 


0.0 


0.003 


0.06 


0.43 


2.12 


1743 


RG0.9c_RP0.24 


0.24 


0.0 


0.003 


0.16 


0.94 


1.41 


261 


RG0.9c_RP0.96 


0.96 


0.0 


0.003 


0.14 


0.77 


1.67 


213 


RG0.9c_RP1.91 


1.91 


0.0 


0.003 


0.10 


0.65 


2.08 


200 


RG0.9c_RP3.82 


3.82 


0.0 


0.005 


0.11 


0.49 


2.79 


325 


RG0.9c_RP5.73 


5.73 


0.0 


0.005 


0.10 


0.36 


3.17 


416 


RG0.9c_RP7.64 


7.64 


0.0 


0.005 


0.12 


0.30 


3.37 


303 



Note. — Final orbital properties and data for the collisions, each labeled in col. (1). Col. (2) shows the 
periastron separation r p , which ranges from nearly head-on to grazing. In col. (3) we find the range of mass 
remaining bound to the point mass representing the subgiant or red giant core: AM C = M 2 — m pt . Similarly, 
col. (4) displays the final mass bound to the NS: AM NS = Mi - 1.4M . Col. (5) gives the mass M 3 of the 
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